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ABSTRACT 


Observations have found black holes spanning ten orders of magnitude in mass across most of cosmic 
history. The Kerr black hole solution is however provisional as its behavior at infinity is incompatible 
with an expanding universe. Black hole models with realistic behavior at infinity predict that the grav- 
itating mass of a black hole can increase with the expansion of the universe independently of accretion 
or mergers, in a manner that depends on the black hole's interior solution. We test this prediction 
by considering the growth of supermassive black holes in elliptical galaxies over 0 < z < 2.5. We 
find evidence for cosmologically coupled mass growth among these black holes, with zero cosmological 
coupling excluded at 99.9896 confidence. The redshift dependence of the mass growth implies that, 
at z S 7, black holes contribute an effectively constant cosmological energy density to Friedmann's 
equations. The continuity equation then requires that black holes contribute cosmologically as vacuum 
energy. We further show that black hole production from the cosmic star formation history gives the 
value of Qa measured by Planck while being consistent with constraints from massive compact halo 
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objects. We thus propose that stellar remnant black holes are the astrophysical origin of dark energy, 
explaining the onset of accelerating expansion at z ~ 0.7. 


Keywords: Supermassive black holes (1663) — Astrophysical black holes (98) — Dark energy (351) 


1. INTRODUCTION 


Astrophysical black holes (BHs), with masses span- 
ning a few to several billion solar masses, are found in 
systems ranging from stellar binaries to supermassive 
BHs in the centers of galaxies. Observations of gravi- 
tational waves from binary BH mergers (Abbott et al. 
2019, 2021) and of supermassive BHs by the Event Hori- 
zon Telescope Collaboration et al. (2019) and Akiyama 
et al. (2022), have shown excellent consistency with the 
Kerr (1963) solution on timescales from milliseconds to 
hours, and spatial scales of up to milliparsecs. BHs are 
thus established as a universal phenomenon, across at 
least ten orders of magnitude in mass. 

Existing models for astrophysical BHs are necessarily 
provisional. They feature singularities, horizons, and 
unrealistic boundary conditions! (e.g. Wiltshire et al. 
2009). Though singularities and horizons are of the- 
oretical interest (e.g. Harlow 2016), the Kerr solution 
reduces to flat spacetime at spatial infinity. This is in- 
compatible with our universe, which is in concordance 
with a perturbed Robertson-Walker (RW) cosmology to 
sub-percent precision (e.g. Aghanim et al. 2020; Dodel- 
son & Schmidt 2020). Thus, regardless of singularities 
and horizons, Kerr is only appropriate for intervals of 
time short compared to the reciprocal expansion rate of 
the universe, and can only be consistently interpreted as 
an approximation to some more general solution. 

Efforts to construct a BH model in General Relativity 
(GR) with realistic RW boundary conditions have been 
ongoing for nearly a century, but have met with lim- 
ited success. Early work by McVittie (1933) generalized 
the Schwarzschild solution to arbitrary RW spacetimes. 
Nolan (1993) constructed a non-singular interior for this 
solution, and progress has been made in understanding 
its horizon/causal structure (e.g. Kaloper et al. 2010; 
Lake & Abdelqader 2011; Faraoni et al. 2012; da Silva 
et al. 2013). Faraoni & Jacques (2007) constructed solu- 
tions featuring dynamical phenomena such as horizons 
that comove with the universe's expansion, evolution of 
interior energy densities and pressures, and time-varying 
mass. These solutions are significant, because they show 
how heuristic application of Birkhoff's theorem in cos- 


mological settings can fail in the presence of strong grav- 
ity? (c.f. Lemaitre 1931; Einstein & Straus 1945; Callan 
et al. 1965; Peebles 1993). Time-varying mass in partic- 
ular has been studied by Guariento et al. (2012) and Ma- 
ciel et al. (2015), but its interpretation remains largely 
unexplored. All of these solutions, however, are incom- 
patible with Kerr on short time-scales because they do 
not spin. A BH solution that satisfies observational con- 
straints at small and large scales simultaneously has yet 
to be found. 

Progress on these problems in GR has become possible 
with advances that resolve a long-standing ambiguity in 
Friedmann’s equations (e.g. Ellis & Stoeger 1987). In 
RW cosmology, the metric is position-independent and 
has no preferred directions in space. In order for Ein- 
stein’s equations to be consistent, the stress-enegy must 
share these properties. Einstein's equations, however, 
give no prescription for converting the actual, position- 
dependent, distribution of stress-energy observed at 
late-times into a position-independent source. Croker & 
Weiner (2019) resolved this averaging ambiguity, show- 
ing how the Einstein-Hilbert action gives the necessary 
relation between the actual distribution of stress-energy 
and the source for the RW model. 

A consequence of this result is that relativistic mate- 
rial, located anywhere, can become cosmologically cou- 
pled to the expansion rate. This has implications for 
singularity-free BH models, such as those with vac- 
uum energy interiors (e.g. Gliner 1966; Dymnikova 1992; 
Chapline et al. 2002; Mazur & Mottola 2004; Lobo 2006; 
Mazur & Mottola 2015; Dymnikova & Galaktionov 2016; 
Posada 2017; Posada & Chirenti 2019; Beltracchi & 
Gondolo 2019). The stress-energy within BHs like these, 
and therefore their gravitating mass, can vary in time 
with the expansion rate. The effect is analogous to cos- 
mological photon redshift, but generalized to timelike 
trajectories. 

The presence or absence of cosmologically coupled 
mass in BHs strongly constrains observationally viable 
BH solutions. In general, the way in which a BH’s mass 
M changes in time depends on the BH model. Cro- 


1 Boundary conditions include both behavior at junctions and 
asymptotic behavior at infinity, if applicable. 


2 It is not sufficient that densities and pressures exterior to some 
region justify a Newtonian treatment. Densities and pressures in- 
terior to that region must also remain non-relativistic (Weinberg 
2008). 
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Figure 1. (Top) Posterior distributions of cosmological coupling strength k, inferred by comparing SMBHs in local elliptical 
galaxies to those in five samples of elliptical galaxies at z > 0.7. (Bottom) Combining these posterior samples with equal 
weighting gives a distribution with k = audies at 9096 confidence. If fit to a Gaussian, the fit has a mean of k — 3.09 with 
a standard deviation of 0.76 (shading). Vertical lines indicate: k — 0 coupling, as expected for traditional BHs like Kerr and 
the decoupled solution by Nolan (1993); and k = 3 coupling, as predicted for vacuum energy interior BHs. The measurement 
disfavors zero coupling at 99.98% confidence and is consistent with BHs possessing vacuum energy interiors, as first suggested 


by Gliner (1966). 


ker et al. (2021) give a parameterization of the effect in 
terms of the RW scale factor a, 


k 
ua) = Mta) (£) a 2 di, (1) 
i 
where a; is the scale factor at which the object becomes 
cosmologically coupled and k > 0 is the cosmological 
coupling strength.” The Nolan (1993) solution can be 
regarded as cosmologically coupled with k — 0 because 
its stress-energy evolves such that the mass remains 
fixed throughout the cosmological expansion. Vacuum 
energy interior solutions with cosmological boundaries 
have been predicted to produce k ~ 3, which is the 
maximum value for causal material with positive energy 
density (Croker & Weiner 2019). 
Cosmologically coupled mass change allows for exper- 
imental distinction between singular and non-singular 


3 This model agrees with the GR prediction for a population of 
identical objects, with homogeneous and non-singular interiors. 
For such objects, k = —3P/p, where P and p are the objects 


interior pressure and energy density, respectively. 


BHs, complementing constraints from short time-scale 
data (e.g. Sakai et al. 2014; Uchikata et al. 2016; Yunes 
et al. 2016; Cardoso et al. 2016; Cardoso & Pani 2017; 
Chirenti 2018; Konoplya et al. 2019; Maggio et al. 2020). 
Observing cosmologically coupled mass change, how- 
ever, is challenging. Between an initial scale factor a; 
and a final one af, mass evolution only becomes appar- 
ent when (ar/a;)" >> 1. For example, an observational 
test of cosmological mass change at z < 3 requires a 
population of BHs whose masses can be tracked across 
at least a Gyr, and in which accretion or merging can 
be independently estimated. 

In this paper, we perform a direct test of BH mass 
growth due to cosmological coupling. A recent study by 
Farrah et al. (2023) compares the BH masses Mgy and 
host galaxy stellar masses M, of ‘red-sequence’ elliptical 
galaxies over 6 — 9 Gyr, from the current epoch back to 
z~ 2.7. The study finds that the BHs increase in mass 
over this time period by a factor of 8 — 20x relative to 
the stellar mass. The growth factor depends on red- 
shift, with a higher factor at higher redshifts. Because 
SMBH growth via accretion is expected to be insignif- 
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icant in red-sequence ellipticals, and because galaxy- 
galaxy mergers should not on average increase SMBH 
mass relative to stellar mass, this preferential increase 
in SMBH mass is challenging to explain via standard 
galaxy assembly pathways (Farrah et al. 2023, 85). We 
here determine if this mass increase is consistent with 
cosmological coupling, and, if so, the constraints on the 
coupling strength k. 


2. METHODS 


We consider five high-redshift samples, and one lo- 
cal sample, of elliptical galaxies given by Farrah et al. 
(2023). For the high-redshift samples we use: two from 
the WISE survey (one at Z — 0.75 measured with the 
H8 line, and one at Z = 0.85 measured with the Mg II 
line), two from the SDSS (one at Z — 0.75 and one at 
Z = 0.85, with Hf and Mg II, respectively) and one from 
the COSMOS field (at Z = 1.6). We then determine the 
value of k needed to align each high redshift sample with 
the local sample in the Mgy — M, plane. If the growth 
in BH mass is due to cosmological coupling alone, re- 
gardless of sample redshift, the same value of k will be 
recovered. 

To compute the posterior distributions in k for each 
combination, we apply the pipeline developed by Farrah 
et al. (2023), which we briefly summarize. Realizations 
of each galaxy sample are drawn from the sample with 
its reported uncertainties. The likelihood function ap- 
plies the expected measurement and selection bias cor- 
rections to the realizations, as appropriate for each sam- 
ple. The de-biased, and so best actual estimate, BH 
mass of each galaxy is then shifted to its mass at z = 0 
according to Equation 1 with some value of k. Using 
the Epps-Singleton test, an entire high-redshift realiza- 
tion is then compared against a realization of the local 
ellipticals, where BH masses are shifted to z — 0 in the 
same way. The result is a probability that can be used 
to reject the hypothesis that the samples are drawn from 
the same distribution in the Mgy — M. plane, i.e. that 
they are cosmologically coupled at this k. 


3. RESULTS & DISCUSSION 


We present posterior distributions in k, for each high- 
redshift to local comparison, in the top row of Figure 
1. The redshift dependence of mass growth translates 
to the same value k ~ 3 across all five comparisons, as 
predicted by growth due to cosmological coupling alone. 
As further verification, we compute k from a compari- 
son between high-redshift WISE and COSMOS samples. 
This comparison requires no BH bias corrections. We 


find a consistent value of k — 06 Combining the 


results from each local comparison gives, 


k =3.111}33 (90% confidence), (2) 
which excludes k = 0 at 99.98% confidence, equivalent 
to > 3.90 observational exclusion of the singular Kerr 
interior solution. 


3.1. Implications 


Our result provides a single-channel explanation for 
the disparity in SMBH masses between local ellipticals 
and their 7-10 Gyr antecedents (Farrah et al. 2023). Fur- 
thermore, the recovered value of k ~ 3 is consistent with 
SMBHs having vacuum energy interiors. Our study thus 
makes the existence argument for a cosmologically re- 
alistic BH solution in GR with a non-singular vacuum 
energy interior. 

Equation (1) implies that a population of k ~ 3 BHs 
will gain mass proportional to a?. Within a RW cos- 
mology, however, all objects dilute in number density 
proportional to a7. When accretion becomes subdom- 
inant to growth by cosmological coupling, this popu- 
lation of BHs will contribute in aggregate as a nearly 
cosmologically constant energy density. From conserva- 
tion of stress-energy, this is only possible if the BHs also 
contribute cosmological pressure equal to the negative of 
their energy density, making k ~ 3 BHs a cosmological 
dark energy species. 


3.2. Qa from the cosmic star formation history 


If k ~ 3 BHs contribute as a cosmological dark en- 
ergy species, a natural question is whether they can 
contribute all of the observed Q4. We test this by as- 
suming that: (1) BHs couple with k = 3, consistent with 
our measured value; (2) BHs are the only source for Qa, 
and (3) BHs are made solely from the deaths of mas- 
sive stars. Under these assumptions, the total BH mass 
from the cosmic history of star formation (and subse- 
quent cosmological mass growth) should be consistent 
with Qa = 0.68. 

In Appendix A we construct a simple model of the 
cosmic star formation rate density (SFRD) that allows 
exploring combinations of stellar production rate, stel- 
lar IMF, and accretion history. Figure 2 displays models 
that produce the Planck measured value of Qa = 0.68 
(Aghanim et al. 2020) with the indicated IMF and Ed- 
dington ratio Ae. Any monotonically decreasing path 
inside the filled region produces f?4 = 0.68 for some ob- 
servationally viable IMF and accretion history, consum- 
ing at most 396 of baryons. This baryon consumption is 
compatible with the results of Macquart et al. (2020), 
who find that Q, at low redshift agrees with Q, inferred 
from the Big Bang to within 50%. 
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Figure 2. Cosmic star formation rate densities (SFRDs) capable of producing the necessary k = 3 BH density to give 
Qa = 0.68 (green, solid). The details of the model are given in Appendix A. The upper bound of the viable region adopts 
a Kroupa (2002) IMF at all redshifts with the least amount of remnant accretion required to produce Qa with a decreasing 
power-law SFRD model (red, solid). The lower bound adopts the top-heavy IMF of Harikane et al. (2022a) at z > 7 (blue, 
solid). Two middle lines show the impact of a top-heavy IMF at z > 7, but no remnant accretion (green, solid); and higher 
accretion, but with a Kroupa IMF (orange, solid). Existing measurements of the SFRD via IR (Rowan-Robinson et al. 2018, 
red, squares), y-ray bursts (Kistler et al. 2009, orange, stars), FIR (Algera et al. 2023, brown, xs), and rest-frame UV via JWST 
(Donnan et al. 2022; Harikane et al. 2022a, purple, dots), (Bouwens et al. 2022, blue, dots) are over-plotted. The UV points can 
vary ^ —1 dex depending upon IMF assumptions and UV luminosity integration bounds. Consistency occurs with consumption 
of < 3% of the baryon fraction Q, after cosmic dawn. The results assume stellar first light at z, = 25 (Harikane et al. 2022b, 


Fig 25) but are typical of the scenario for 15 < z, < 35. Also indicated are the redshifts probed by 21cm experiments. 


It follows from Equation (1) that cosmological cou- 
pling in BHs with & = 3 will produce a BH population 
with masses > 10?Mo. If these BHs are distributed in 
galactic halos, they will form a population of MAssive 
Compact Halo Objects (MACHOs). In Appendix B, we 
consider the consistency of SFRDs in Figure 2 with MA- 
CHO constraints from wide halo binaries, microlensing 
of objects in the Large Magellanic Cloud (LMC), and 
the existence of ultra-faint dwarfs (UFDs). We conclude 
that non-singular k = 3 BHs are in harmony with MA- 
CHO constraints while producing Q4 = 0.68, driving 
late-time accelerating expansion. 


4. FUTURE TESTS 


Further tests of cosmological coupling in BHs are es- 
sential to confirm or refute our proposal. We list some 
examples below, highlighting possible tensions. 


4.1. Signatures in the cosmic microwave background 


A population of k ~ 3 BHs are a dark energy species. 
Thus, their distribution in space need not trace baryons 
or dark matter at all times. For example, Croker et al. 


(2020b) study the spatial distribution of k ~ 3 BHs 
with anisotropic stress (e.g. Cattoen et al. 2005) at first 
order in cosmological perturbation theory. Anisotropic 
stress within individual BHs leads to an effective fluid 
at first-order that resists clustering and can even drive 
the spatial distribution to uniformity.4 Resistance to 
clustering is computed from the non-singular BH model, 
of which there is currently none preferred (84.6). Some 
amount of anisotropic stress is, however, necessary to 
satisfy constraints on the galaxy two-point correlation 
function (Croker et al. 2020b). 

Anisotropic stress from cosmologically coupled BHs 
at z, = 25 would maximally impact the CMB via the 
Integrated Sachs-Wolfe effect at 4 = 5, with lesser im- 
pacts at £ = 5 (e.g. Koivisto & Mota 2006; Dodelson 
& Schmidt 2020). The low £ < 30 of the CMB TT 
anisotropy spectrum are anomalous in several respects 
(e.g. Schwarz et al. 2016) and with amplitude in excess 
of cosmic variance at £ — 5 (e.g. Planck Collaboration 


4 For mathematical details, see Croker et al. (2022) 
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et al. 2020). Any imprint of BH production at cosmic 
dawn on the CMB low ¢’s provides a precision test of 
non-singular BH models, and could enable independent 
constraint of the high-z SFRD. 


4.2. Strong lensing of y-ray bursts 


'The prevalence of strong gravitationally lensed GRBs 
has been used to estimate the comoving density of In- 
termediate Mass Black Holes (IMBHs) (Paynter et al. 
2021, see also e.g. Wang et al. 2021; Yang et al. 2021; 
Chen et al. 2022; Lin et al. 2022; Kalantari et al. 2022; 
Liao et al. 2022). GRBs probe IMBHs at scales much 
greater than the non-linearity scale of ~ 6 Mpc, comple- 
menting ~ 50 kpc MACHO constraints (Appendix B). 
At these larger scales, the prediction of uniform spatial 
distribution from anisotropic stress (§4.1) is directly ap- 
plicable. 

Xiao & Schaefer (2011) report GRBs with spectro- 
scopic redshifts and find that 90% of GRBs have red- 
shifts z > 0.25, with 50% occurring at z > 1.5. Across 
these distances, k ~ 3 cosmological coupling of BHs 
strongly impacts the comoving mass density of lenses 
along the line-of-sight, reducing it by factors of at least 
2 — 16x. Existing analyses do not incorporate this 
impact of cosmological coupling on the optical depth. 
Thus, a direct comparison of measured Opypg against 
the simulated BH populations in Appendix B is not yet 
possible. A promising future test is to incorporate cos- 
mological coupling into the optical depth calculation. 
In general, a spectroscopic redshift is not available for 
candidate lensed GRBs, significantly weakening mea- 
surement of Opugg (e.g. Paynter et al. 2021). Consis- 
tency with our population predictions, when evaluated 
at the lens redshift, could help to constrain such mea- 
surements. 


4.3. Stellar mass BH-BH merger rates 


The population properties of observable binary BH 
mergers probe cosmological coupling because their in- 
spiral time can be a significant fraction of the age of the 
universe. If k > 0, then the component masses of bi- 
nary BH systems observed by gravitational-wave detec- 
tors are not representative of the birth masses. This im- 
pacts the interpretation of merging populations of BHs 
(The LIGO Scientific Collaboration et al. 2021) and the 
stochastic background (Arzoumanian et al. 2020; Chris- 
tensen 2019). Coupled mass growth further leads to 
accelerated orbital decay, though the rate of this decay 
may be model dependent (c.f. Croker et al. 2020a; Had- 
jidemetriou 1963). This aspect affects both the observed 
mass spectrum and rate of binary BH mergers. For ex- 
ample, in the absence of cosmological coupling, stellar- 
mass binary BH systems with semi-major axis = 0.3 AU 


can only merge within a Hubble time when highly ec- 
centric. Conversely, with k > 0.2, systems with initial 
semi-major axes 0.1 < R < 10* AU can merge within 
less than a Hubble time. This can lead to significant 
increases in merger rate. Competing with this effect, 
however, is that mergers can happen so quickly after 
remnant formation that they merge at a redshift beyond 
the detection limit of current generation observatories. 


4.4. Direct measurement from orbital period decay 


Direct measurement of altered orbital decay rate in 
binary BH systems with proposed space-based observa- 
tories like LISA (Armano et al. 2018), DECIGO (Sato 
et al. 2017), and Taiji (Gong et al. 2021) is likely not 
possible because these observatories can only track or- 
bital decay in the final few months. Decade-scale elec- 
tromagnetic observations of a pulsar-BH orbit, however, 
would provide sufficient precision to measure cosmolog- 
ical coupling directly (e.g. Croker & Weiner 2019; Weis- 
berg et al. 2010). 


4.5. Stellar mass BHs and stellar evolution 


Mass shifts consistent with cosmological coupling have 
been proposed to exist in the merging binary BH pop- 
ulation, explaining both the observed BH mass spec- 
trum and the existence of BHs in the pair-instability 
mass gap, though with a smaller coupling strength of 
k ~ 0.5 (Croker et al. 2021). A coupling of k = 3 
and adopting contemporary stellar population synthe- 
sis estimates can lead to an overabundance of BHs with 
masses > 120Mo. While uncertainties in binary BH for- 
mation channels (Mandel & Farmer 2022; Zevin et al. 
2021), massive binary star evolutionary physical pro- 
cesses (Broekgaarden et al. 2022), nuclear reaction rates 
(Farmer et al. 2020), supernova core collapse physics 
(Patton & Sukhbold 2020), and SFRD and metallic- 
ity evolution (Chruslitiska 2022) leave population model 
flexibility, there are known young BHs within X-ray bi- 
naries with mass ~ 20Mo (e.g. Miller-Jones et al. 2021). 
If this BH mass is typical of young stellar remnants at 
z S 5, then the distribution of remnant binary semi- 
major axes and eccentricities becomes constrained so as 
not to produce overly massive BH-BH mergers. An im- 
portant test for k = 3 BHs is whether such constraints 
are plausible. 


4.6. Theoretical considerations 


As described in 81, there are known exact solutions 
with each of the following properties: strong spin, ar- 
bitrary RW asymptotics, dynamical mass, and interior 
vacuum energy equation of state. Our result implies the 
existence, within GR, of an exact solution with all of 


COSMOLOGICAL COUPLING OF BLACK HOLES: FIRST OBSERVATION AND IMPLICATIONS 7 


these properties. Currently, there is no known solution 
that possesses all four, though there are known solu- 
tions with various combinations of two (e.g. Guariento 
et al. 2012; Dymnikova & Galaktionov 2016). Finding 
solutions that feature all four properties is an important 
theoretical step forward. 

It is also interesting that a link between BHs and late- 
time accelerating expansion has been independently sug- 
gested within frameworks distinct from GR. Afshordi 
(2008) and Prescod-Weinstein et al. (2009) adopt a grav- 
itational æther framework in which an effective cosmo- 
logical constant emerges from quantum gravity effects 
at BH horizons. Notably, this yields a reduction of 
the quantum field theory tuning required, from over one 
hundred decimal places to only two. Their scenario is 
however distinct from ours; it has a different theoretical 
basis and does not, to our knowledge, predict that BHs 
gain mass as the scale factor increases. 


4.7. Validation of preferential growth of SMBHs 


Further validation of the measured preferential growth 
of SMBHs by Farrah et al. (2023) is an essential test of 
our proposal. This can be done by: deepening under- 
standing of the relevant biases, improving morphologi- 
cal determinations of the high redshift samples, testing 
combinations of traditional assembly pathways in simu- 
lations, and improving the accuracy of SMBH and stellar 
mass measures. An important future test is to improve 
the statistics and reliability of the high redshift sample. 
Doing so requires assembling a sample of > 10° AGN 
in elliptical hosts with low SFRs and reddenings over 
0.7 < z < 2.5. The sample should have reliable measures 
of host stellar mass and consistent measures of SMBH 
mass with a subset that includes multi-epoch reverber- 
ation mapped measures. Such a study would enable 
narrower redshift intervals to be used at z > 0.7, giving 
finer discrimination between cosmological coupling and 
other processes. 


4.8. Quasars at z > 6 and the SMBH mass function 


The existence of SMBHs at z = 6 (e.g. Trakhtenbrot 
et al. 2015) with masses > 10°Mo is challenging to ex- 
plain via accretion and direct collapse models (e.g. In- 
ayoshi et al. 2020; Volonteri et al. 2021). Cosmological 
coupling with k = 3, starting at z = 25, provides a 
mass increase of a factor of 51 by z = 6 (Equation 1). 
This would ease tensions between BH growth models 
and observations of z > 6 quasars, but it has not been 
shown that cosmological coupling is required to do so. 
BH masses must also increase between z ~ 6 and z = 0 
by a factor of 343. Comparison of the BH mass func- 
tion in quasars at z ~ 6 with the BH mass function at 


z = 0 must be compatible with this minimum increase. 
Generalization of the Soltan (1982) argument to k = 3 
coupling is a first step, though comparison of the high- 
end mass cutoff is not sufficient, because SMBHs may 
cease to accrete luminously above ~ 101! Mc (e.g. King 
2015; Inayoshi & Haiman 2016; Carr et al. 2021). 


5. CONCLUSIONS 


Realistic astrophysical BH models must become cos- 
mological at large distance from the BH. Non-singular 
cosmological BH models can couple to the expansion of 
the universe, gaining mass proportional to the scale fac- 
tor raised to some power k. A recent study of supermas- 
sive BHs within elliptical galaxies across ~ 7 Gyr finds 
redshift-dependent 8 — 20x preferential BH growth, rel- 
ative to galaxy stellar mass. We show that this growth 
excludes decoupled (k — 0) BH models at 99.9896 con- 
fidence. Our measured value of k = 3.11*119 at 9096 
confidence is consistent with vacuum energy interior BH 
models that have been studied for over half a century. 
Cosmological conservation of stress-energy implies that 
k — 3 BHs contribute as a dark energy species. We show 
that k — 3 stellar remnant BHs produce the measured 
value of Q4 within a wide range of observationally viable 
cosmic star formation histories, stellar IMFs, and rem- 
nant accretion. They remain consistent with constraints 
on halo compact objects and they naturally explain the 
"coincidence problem," because dark energy domination 
can only occur after cosmic dawn. Taken together, we 
propose that stellar remnant k — 3 BHs are the astro- 
physical origin for the late-time accelerating expansion 
of the universe. 
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APPENDIX 


A. SFRD BOUNDING MODEL 


Here we construct a simple model to establish whether 
plausible SFRDs lead to Qa œ 0.7, if massive stellar 
collapse produces k = 3 BHs. For the overall form of the 
SFRD, we adopt the model of Madau & Fragos (2017) 
at z < 4. Between z = 4 and cosmic dawn at z = 25 we 
adopt a power-law for the SFRD with exponent o « 0, 
matched continuously to the z < 4 SFRD. In units of 
Mo: Mpc? eyri, 


(14-2)?:6 


(2/4)? p, (o) 4«z«25. (Al) 
0 z 25 


p«(z,a) = 


where p, is the time rate of production of stellar mass. 
The choice of varying the power law slope at z > 4 is mo- 
tivated by the disparities in observed p, over 4 < z < 10. 
Estimates from infrared (Rowan-Robinson et al. 2018) 
and y-ray bursts (Kistler et al. 2009) can be 1-2dex 
higher than estimates from UV-dropouts (Donnan et al. 
2022; Harikane et al. 2022a; Bouwens et al. 2022). Vary- 
ing a at z > 4 allows us to encompass SFRDs consistent 
with all extant data. 

To account for variation in the stellar IMF with red- 
shift, we divide the SFRD into a Population III (Pop 
III) epoch at z > 7 (e.g. Inoue et al. (2014) but c.f. 
Liu & Bromm (2020)) and a standard epoch, character- 


ized by a Kroupa IMF at z < 7. We consider scenarios 
where the Pop III epoch is characterized by either a 
Kroupa IMF or the top-heavy Harikane et al. (2022a) 
IMF œ M-??? for 50 < M < 500, and zero elsewhere. 
Finally, we use the redshift-dependent, mean-metallicity 
model of Madau & Fragos (2017), truncated at Zo, to 
enable use of standard metallicity-dependent zero-age 
main sequence (ZAMS) mass to remnant mass, delayed, 
core-collapse supernovae prescriptions by Fryer et al. 
(2012). 

For simplicity, we apply post-remnant formation ac- 
cretion, with a fixed Eddington ratio Ae, to all Pop III 
BHs over a duration te, 


M(t; + te)  M(t;) exp (Aete/€), (A2) 


brief enough that the impact on instantaneous rate from 
cosmological coupling of accreted mass can be neglected. 
Because te, Ac, and e are degenerate in Equation (A2), 
there is effectively only one parameter. Farrah et al. 
(2022) use simulation studies of galaxy mergers to esti- 
mate t, = 44+22 Myr at Ae > 1 from observations of 21 
infrared-luminous galaxy mergers at z < 0.3. This value 
is consistent with ~ 10 Myr estimates by Madau et al. 
(2014) for high z BH seeds, so we adopt t, = 20 Myr. 
The uncertainty in te is large, e.g. Safarzadeh & Haiman 
(2020) find Pop III stellar BH accretion across timescales 
^ 0.5 Myr. Tortosa et al. (2022) measure Ae ~ 472 for 
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a Seyfert 1 galaxy, while simulations by Inayoshi et al. 
(2016) find stable accretion onto isolated Pop III BHs 
at Ae ~ 5000. As we apply this accretion to all Pop III 
BHs in our model, we consider models within a more 
conservative Ae < 30. Given the large range of plausible 
te and Ae, we fix the efficiency e = 0.16, roughly 30% of 
the Kerr BH limit (e.g. Bardeen 1970). 

'To compute viable SFRDs, we use a Monte Carlo ap- 
proach. Because high-mass remnants formed at high 
redshift can easily dominate Q,, care must be taken to 
sufficiently sample these tails of the distributions. We 
divide the domains for IMF and SFRD distributions into 
distinct windows such that, in these windows, the dis- 
tributions are guaranteed to give Poisson errors < 1% 
in the lowest probability bin. The total collection of 
draws is eventually re-weighted by relative areas under 
the windowed IMFs and Equation A1. In each window, 
we draw 10° ZAMS stellar masses from the appropriate 
IMFs and 10° redshift-dependent metallicities. Given 
the birth masses and redshifts, we approximate a red- 
shift of stellar death using für, ~ 10!°(M/Mo)~7 yr, 
where 7 — 2.5 (e.g. Harwit 2006). We discard draws that 
live beyond z — 0, and then determine collapsed rem- 
nant masses. The ZAMS mass to remnant mass model 
we adopt assumes that all stars drawn are single stars. 
Consistent with Kalogera & Baym (1996), we regard any 
remnants with mass > 2./7M as BHs. To convert the 
resulting population of BHs into a cosmological density, 
we regard this drawn population as residing within one 
co-moving Mpc?. We first reweight all draws by the 
probability of having come from their respective win- 
dows. We then convert from this reweighted mass to a 
predicted mass by rescaling p.(z = 0) by the total mass 
in long-lived stars, neutron stars, and initial BH rem- 
nant mass, as described by Madau & Fragos (2017). We 
compute „(z = 0) by integrating Equation A1 in time, 
and scaling by an IMF-appropriate gas return fraction 
1— R (e.g. R = 0.39 given by Madau & Fragos (2017) 
for Kroupa, and 1 for Harikane). We then divide the 
summed predicted mass in BHs at z — 0 by the critical 
density today 3H} /87G to get Qa. 


B. MASSIVE COMPACT HALO OBJECT 
CONSTRAINTS 


Here we establish that a k ~ 3 BH population suf- 
ficient to produce Q; is consistent with constraints on 


massive compact halo objects (MACHOs). In Figure 3, 
we display the contribution to halo density from k = 3 
BH masses at z = 0 computed in the explicit SFRD 
models shown in Figure 2. Consistent with the discus- 
sion in 84.1, we have assumed a uniformly dispersed pop- 
ulation in the computation of Figure 3. We adopt a uni- 
form Milky Way (MW) halo density equal to the median 
DM density at the Solar System 8.8+0.5 x 107? Mo /pc? 
as measured by Cautun et al. (2020). This is conserva- 
tive relative to a mean halo density 1.7 x 107? Mo /pc?, 
inferred from 1.37 x 10'! Mc within the MW at < 20 kpc 
as measured by Posti & Helmi (2019). 

We display limits on MACHOs from microlensing (e.g. 
Blaineau et al. 2022), wide halo binary disruption (e.g. 
Tyler et al. 2022; Monroy-Rodríguez & Allen 2014), and 
UFD disruption (e.g. Brandt 2016). We approximate 
the UFD constraint from Brandt (2016) as c 1/m fol- 
lowing Binney & Tremaine (2011, Eqn. (7.104)) and 
convert from the UFD halo density assumed by Brandt 
(2016, Figure 4) to our MW halo density scale by a mul- 
tiplicative factor 0.3/0.0088. 

All BH density contributions lie below the microlens- 
ing and wide binary limits by several orders of magni- 
tude. High accretion models with a Kroupa IMF show 
some tension with UFD constraints, which begin to have 
discriminatory power in our scenario. Comparing the 
two top-heavy IMF models, the effect of accretion is to 
redistribute mass density from the IMBH range into the 
SMBH range. The top-heavy IMF with no accretion 
model shows how a top-heavy IMF acts to populate the 
IMBH region. 

The aforementioned constraints do not account for co- 
moving BH mass density decrease with increasing red- 
shift (Equation (1)) for k > 0. The constraints as pre- 
sented are thus overly stringent because the probabil- 
ity of compact object interactions is reduced at earlier 
times. The UFD galaxy and wide halo binary con- 
straints are most impacted, as they consider the stabil- 
ity of astrophysical phenomena over Gyr timescales. A 
thorough analysis that incorporates decreasing comov- 
ing density of compact objects further back in time is 
the subject of future work. 
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